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Abstract 

Non-linear, dynamic, finite element analysis is used in various engineering disciplines to evaluate high-speed, 
dynamic impact and vibration events. Some of these applications require connecting rotating to stationary 
components. For example, bird impacts on rotating aircraft engine fan blades are a common analysis performed 
using this type of analysis tool. Traditionally, rotating machines utilize some type of bearing to allow rotation in 
one degree of freedom while offering constraints in the other degrees of freedom. Most times, bearings are modeled 
simply as linear springs with rotation. This is a simplification that is not necessarily accurate under the conditions 
of high-velocity, high-energy, dynamic events such as impact problems. For this reason, it is desirable to utilize a 
more realistic non-linear force-deflection characteristic of real bearings to model the interaction between rotating 
and non-rotating components during dynamic events. The present work describes a rolling element bearing model 
developed for use in non-linear, dynamic finite element analysis. This rolling element bearing model has been 
implemented in LS-DYNA as a new element, ^ELEMENT BEARING. 

Introduction 

While there are many different types of bearings, all perform the same basic purpose of allowing 
rotation in one degree of freedom, while constraining motion in the others. Most bearings can 
also be separated into two major classes: fluid film (or journal) bearings and rolling element 
bearings. Fluid film bearings are defined by a film of fluid (liquid or gas) separating the rotating 
and non-rotating components. Fluid film bearings are common in applications such as large 
industrial steam turbines, industrial compressors, pumps, blowers, and similar types of machines. 
Rolling element bearings describe a class that includes any bearing with rolling elements 
(cylindrical, spherical, and conical rollers, balls, and needles for example) separating the rotating 
and non-rotating components. Rolling element bearings are common in applications such as 
aeronautic and automotive engines, wheel bearings, gearboxes, and the like. Both classes of 
bearings have another common feature; they exhibit a non-linear relationship between the 
relative deflection of the rotating and non-rotating parts, and the load applied. One could say 
they behave like non-linear springs. Often, but not always, bearings are stiff compared to the 
structures to which they are connected, or loads are relatively light and therefore, deflections are 
typically small. In these cases, it can be appropriate to model them with linear load-deflection 
behavior, i.e. like linear springs. This is a very co mm on practice in rotordynamics and finite 
elements analysis. However, in the cases where the bearings are soft compared to the 
surrounding structure, or loads are very large, the non-linearity of the bearings can be very 
important. The motivation for the current work is to more appropriately model bearings, 
specifically rolling element bearings, in finite elements models for cases where the non-linear 
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characteristics are important. Two common examples are aero engine blade-out events and bird 
ingestion events where the loads on the bearings can be very high. It is anticipated that the 
resulting non-linear bearing load-deflection model will be useful for other types of analyses as 
well. 


Contact Mechanics for Bearings 

In general, rolling element bearings consist of four basic components: rolling elements 
(commonly spherical balls or cylindrical rollers), an inner ring (or race), an outer ring (or race), 
and a cage. The function of a bearing is to support radial, axial, and moment loads while 
providing very low resistance about their axis of rotation. As loads are applied, the inner ring 
and outer ring move relative to each other, causing the rings to compress one or more of the 
rolling elements. The races and the rolling elements deform due to the contact pressure at each 
contact interface. The global load-deflection relationship of the bearing is the composite of the 
individual contact interface load-deflection characteristics. 

For the work he performed in the late 1800’s, Hertz (1) is credited with being the first to study 
what is now referred to as Contact Mechanics - the study of deformable bodies in contact with 
one another. The stress that develops in the small contact area of two elastic surfaces has come 
to be known as Hertzian Stress. Others have extended Hertz’s work to estimate not only the 
contact stress, but also the load-deflection relationships in more complicated geometrical 
contacts. As such, representations of the load-deflection relationship in rolling element bearings 
are often based on Hertzian stress calculations. 


The treatment of contact mechanics in rolling element bearings is not new; in fact, the 
fundamentals are described in numerous works on the subject of bearings and contact mechanics, 
some examples of which appear in (2-5). The details of the derivations are left to those texts for 
the interested reader. However, a short summary is given here to describe the basic concepts as 
they are applied to the current bearing model. (The treatment can be extended as shown in (6,7) 
to estimate the stiffness of a ball bearing under a given set of load conditions. This was not done 
for the current model, but this approach can be used to linearize the model for small 
perturbations analysis.) Following the treatment by Harris and Kotzalis (3), which is based on 
Hertz’s work and is similar to the other examples, one can generate relationships for the 
deflection in a contact between two curved bodies, and the normal force applied to push them 
together. In general, for two curved bodies in contact, the curvature sum and curvature 
difference are defined respectively as: 


Zp = 





f(p) 


Ep 


( 1 ) 

( 2 ) 


where, in the case of a spherical ball bearing, r/i = 172 = the ball radius, r//i, r m are the radii of 
curvature of the race groove in the circumferential and axial directions, respectively. 

Consider the contact of the ball and inner raceway in an angular-contact ball bearing shown in 
Fig. 1 . Let body I and II correspond to the ball and the inner race, respectively, and let principal 
direction 1 coincide with the radial plane and principal direction 2 lie transverse to the direction 
of rolling. Using these identities, 
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r ll - r l2 ~ 2 

(3) 

r ni =-(^-D) 

111 2 \cdb a: / 

(4) 

and 


r lIZ = fi® 

(5) 

where / is a factor relating the raceway groove depth to the ball diameter and d m 
bearing pitch diameter. For the outer race contact, r u and r [2 are the same as 
contact, but 

is the nominal 
the inner race 

r ni = +o) 

JJ1 2 ycd s or / 

(6) 

and 


r-j 

II 

1 

(7) 

The initial contact angle, ao, is the contact angle (as shown in Fig. 1) between the balls and the 
races when the bearing is not loaded. In addition, when not loaded, there may be a clearance, Pd, 


between the balls and races. 



Figure 1. Cross section of a typical angular contact ball bearing, adapted from Harris and 
Kotzalis (3). 

One of the primary assumptions made by Hertz, and thus for the method of Harris and Kotzalis, 
is that the shape of the contact between the two curved bodies forms an ellipse. With that 
assumption, the Hertzian formulation gives the curvature difference, Eq. 8, as a function of the 
semi-major (a) and semi-minor ( b ) axes of the contact ellipse, and the two complete elliptical 
integrals of the first and second kind: 
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where 


f0>) = 


(k 2 +i]E-2F 

(k z -1)1 



Li 


K — — 
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( 8 ) 

(9) 

( 10 ) 


( 11 ) 


From here, Harris and Kotzalis show that a relationship can be developed between the applied 
load and the relative approach of the two bodies as: 


[ JS- + (i-yg) \i 2/ ^p 

w V2 k 2 E/ LzSpV ^ )\ 2 


( 12 ) 


where <5= the relative approach of the two bodies in contact (i.e. deflection), Q = applied load, 
V x = Poisson’s ratio of material x, = Young’s modulus of material x. Equation 12 can be 
rearranged into the form: 

Q=KS 3 /z (13) 


where 



(14) 


Equation 14 represents the load-deflection relationship for a single contact (i.e. ball to inner race 
contact or ball to outer race contact). Various methods exist for solution of the elliptical 
integrals (E and F), ranging from lookup tables to numerical solutions. One useful 
approximation is a least squares linear regression developed by Brewe and Hamrock (8) which 
results in errors on the order of less than 5%. 


Now, to determine the composite load-deflection relationship for an entire bearing, Eq. 14 must 
be applied at each contact (i.e. ball to inner race contact and ball to outer race contact) with the 
appropriate load for that contact, and summed over all the rolling elements to give the total load- 
deflection behavior for the bearing. 

To determine the load-deflection relationship for each rolling element, consider a single rolling 
element in a bearing and its contact between the inner race and outer race as shown in Fig. 2. 
One must apply Eq. 14 to the inner race contact and the outer race contact, such that the total 
deflection for one rolling element is: 

d i nn . BT ~F (15) 

Where for the j th rolling element, 5}is the total deflection, and S inngr , S DUtgr are the inner race 
and outer race contact deflections, respectively. 
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Substituting Eq. 14 for each race into Eq. 15 and rearranging gives the overall load-deflection 
relationship for the j th rolling element: 

Q, = ( 16 ) 


where 



Figure 2. Kinematic model of a ball bearing depicting the inner race degrees-of-freedom and 
bearing forces. 


Ball Bearing Load-Deflection Relationship 

Because there are multiple rolling elements in a bearing spaced around the circumference, any 
applied load is shared unevenly by the rolling elements. It can be shown that the total bearing 
forces and moments in response to relative displacement between the inner and outer race are 
found by summing the vector components of the compressive loads and moments from each ball 
contact, giving 


and 



{ COS Ct ; - COS Ipj ) 

cosffjSin^- l 
sin 


(18) 
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Note that the races and the balls have uniform geometry, so the factor is taken out of the 
summation. Also, M z is identically zero in Eq. 19, which indicates that this ideal bearing model 
allows for free axial rotation of the inner ring. 

Implementation in LS-DYNA 

The load-deflection relationship established above for a spherical ball bearing is used to develop 
a non-linear bearing model in LS-DYNA for use in the simulation of dynamic rotating systems. 
At each time step, the model uses the current relative position of two nodes (representing the 
centers of the inner and outer races) and calculates a force to generate a new relative position for 
the next time step, and so on. In this manner, the bearing is represented as a sequence of 
deflections, and corresponding forces and moments, in the global LS-DYNA model for time 
transient solutions. The load-deflection behavior is non-linear, and should be more realistic for a 
bearing subjected to large forces than the more traditional linear spring model. For cases where 
bearing forces are relatively small, the linear spring model may be preferred to minimize 
computational times. 

Implementation of the non-linear bearing load-deflection model in LS-DYNA is accomplished 
through the *ELEMENT_BEARING card. The bearing is represented by two nodes, one 
representing the center of the inner race (usually attached to the rotating part), and the other 
representing the center of the outer race (usually attached to the non-rotating part). Typically, 
these two nodes will be coincident when no load is applied (or when only a preload is applied). 
The user must supply the geometry information in order to calculate a bearing’s load-deflection 
behavior. There are five mandatory cards in *ELEMENT_BEARING. Card 1 defines the 
bearing ID, type of bearing (currently limited to ball bearing), the two nodes, two coordinate IDs 
to define the axes of the rotating and non-rotating nodes, and the number of rolling elements. 
Card 2 defines the Young’s modulus and Poisson’s ratio for both the rolling elements and the 
races, and a flag for yield stress limit. (When the yield stress limit is exceeded a message is 
written to the d3hsp file, and there is no change in the element’s behavior.) Cards 3 and 4 are 
used to define the geometry of the bearing including: the rolling element diameter, bore 
diameter, outside diameter, pitch diameter, initial contact angle, inner and outer race groove 
radius to ball diameter ratios, and the bearing clearance. Lastly, card 5 is used to define any 
preload conditions. There is a preload flag to define the type of preload, and five variables to set 
either x, y, z displacements and x and y rotations, or x, y, z forces and x and y moments. Refer 
to Fig. 1 for a description of the geometry variables. 

Preload is defined as either a set of forces and moments, or a set of displacements and angles that 
mimics the preload in a real system. For example, a typical method of providing preload to 
bearings in a rotor system is to fix the inner race on the rotating shaft by press fit or other 
mechanical means, while the outer race is constrained radially in the non-rotating structure, but is 
free to slide axially. Springs are then used to apply an axial preload force to the outer race to 
prevent the bearings from becoming unloaded and allowing the rolling elements to skid or slide 
at the contact surface with the races, or becoming overloaded and negatively affecting contact 
stress. Skidding is detrimental to bearing life, causing local heating, surface damage, galling, 
increased friction, etc. Therefore, preload is almost always used to apply some nominal contact 
forces. Preload also affects the load-deflection relationship, so it is important to include it in the 
bearing model. 
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Because the bearing load-deflection relationship is non-linear, the stiffness is different depending 
on the load. Therefore, preload effectively changes the stiffness of a given bearing, such that a 
preloaded bearing will behave differently than a non-preloaded bearing. To account for this in 
the DYNA bearing model, the force generated in the bearing is defined as: 

F brs =T(S a +S nodg -)-F prg (20) 

Where S nodg is the instantaneous displacement between the two bearing nodes. F brg is the force 
between the two nodes due to displacement S nodg . 

To account for the modified stiffness from the preload, one must calculate the force due to the 
preload displacement (u 0 ) combined with the instantaneous displacement between the nodes 
(S nodg ), and then subtract the preload force (F pre ) as shown in Eq. 20. Notice, because of the 
non-linearity, this is not equivalent to calculating ^brg from the load-deflection relationship 
directly using S nodg as the displacement. In *ELEMENT_BEARING, Eqs. 18 and 19 are 
implemented to calculate the instantaneous force and moment. If the user supplies preload 
forces, iteration is required to determine iS Q from the load-deflection curve, and a Newton- 
Raphson scheme is used. Alternately, if preload displacements are supplied, F prg is determined 
directly from the load-deflection relationship. In either case, this calculation is only required 
once at the beginning of the simulation, as these values are constant for a given analysis. The 
preload displacement and preload force are considered internal quantities, and are not output as 
any nodal displacements or forces, they are simply required to ensure the instantaneous forces 
and displacements are calculated at the correct location on the load-deflection curve. A side 
benefit of this method is the elimination of the zero force-zero stiffness numerical issues that are 
possible at zero displacement on the load-deflection curve. Users should take note that if they 
choose not to provide preload in their analysis, they may encounter numerical instability if their 
solution involves passing through zero displacement. 

Results and Discussion 

As a check of the model, a standalone computer program utilizing the same model as developed 
here has been written in Visual Basic and implemented within Microsoft® Excel to solve Eqs. 18 
and 19 for a bearing case with documented results useful for comparison. For comparison 
purposes with the published results, the model was extended as mentioned above to calculate 
stiffness from the load-deflection relationship. 

To verify the results from the bearing model, the radial and axial ball bearing stiffness 
coefficients are plotted as a function of radial force and compared to analytical expressions given 
by Gargiulo (6). Fig. 3 shows the results for a bearing whose design parameters are depicted in 
Table 1. The radial stiffness represents the case when the initial contact angle is ao = 0, and the 
axial stiffness represents the case when ao = 90°. The differences between the present model and 
the published results for the radial and axial stiffness coefficients are all less than 3% and 11%, 
respectively, and the trends match qualitatively. 


Table 1 . Design Parameters for Bearing Model Comparison Case 


Elastic modulus of ball 

Eball 

207 

GPa 
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Poisson ratio of ball 

Vball 

0.30 


Elastic modulus of race 

Erace 

207 

GPa 

Poisson ratio of race 

Vrace 

0.30 


Diameter of ball 

D 

19.1 

mm 

Bore diameter 

D, 

100 

mm 

Outer bearing diameter 

D 0 

180 

mm 

Nominal pitch diameter 

d m 

140 

mm 

Initial contact angle 

ao 

40 

deg 

inner raceway radius to ball diameter ratio 

f 

0.54 


outer raceway radius to ball diameter ratio 

fo 

0.52 


Number of balls 

z 

18 


Clearance 

Pd 

0.1 

|um 
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Figure 3. Comparison of the direct radial and axial stiffness coefficients with analytical 
expressions found in literature. 


Conclusions 

A nonlinear formulation has been developed relating the bearing forces and moments to the 
motion of the inner ring in ball bearings. The force-displacement relationships are based on a 
kinematic model of the inner ring with five degrees of freedom. The model accounts for the 
elastic deformation of the balls and races near the contact points and neglects the effects of fluid 
film lubrication and inertia. The model has partially been validated by comparison to closed- 
form expressions found in literature. 

The load-deflection model developed and described herein has been implemented in LS-DYNA 
as *ELEMENT_BEARING and will be available in the next general release of LS-DYNA. The 
model is currently only developed for ball bearings, but future work is planned to extend the 
model to roller bearings and potentially other bearing geometries as needed. 
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